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PREFACE 

The  work  performed  by  the  Communications  group  of  the  University 
of  Alabama  in  Huntsville  for  the  U.  S.  Army  Missile  Command  (MICOM)  under 
contract  DAAH01-71-C-1181,  during  the  period  June  1,  1974  to  July  31,  1975 
consists  of  two  tasks.  The  first  task  is  documented  in  UAH  Research  Report 
No.  176  entitled,  "A  Stvdy  of  the  G-H-K  Tracking  Filter"  and  the  second 
task  will  be  documented  in  a UAH  Research  Report  entitled,  "Quantization 
Effects  in  Digital  MTI  Filters".  This  effort  will  be  denoted  as  fourth- 
year  effort  because  it  is  an  extension  of  a contract  on  EAR  simulation 
studies  started  on  May  18,  1971.  The  goals  of  these  four  consecutive  efforts 
are  reviewed  briefly,  to  place  them  in  the  proper  perspective. 

The  overall  objective  of  the  first  year  effort  was  to  develop, 
implement  and  exercise  a detailed  digital  computer  simulation  of  the  experi- 
mental Array  Radar  (EAR)  system  in  order  to  support  and  complement  the 
experimental  work  being  performed  by  MICOM  personnel.  The  specific  objec- 
tives of  the  study  were  to:  (1)  duplicate  in  a software  form  the  operation 

of  system  hardware  and  its  performance  as  measured  in  the  test  bed,  (2)  to 
simulate  proposed  changes  of  the  hardware  and  predict  performance  improve- 
ment, and  (3)  to  experiment  with  alternate  digital  and  data  processing  tech- 
niques. 

The  objectives  of  the  second  year  effort  were  to:  (1)  refine  and 

update  the  baseline  EAR  simulation  program  and  develop  subprograms  to  study 
critical  components  of  the  system,  (2)  expand  baseline  EAR  simulation  pro- 
grams to  include  multiple  target  handling,  (3)  investigate  via  simulation  a 
constant-false-alarm-rate  (CFAR)  processor,  and  (4)  investigate  the  feasi- 
bility of  glint  error  reduction  through  the  use  of  frequency  agility,  pulse 
compression  or  both. 


ii 


The  objectives  of  the  third-year  effort  were  (1)  A/D  converter 
evaluation  by  digital  computer  analysis  and  (2)  an  expansion  of  the  scope 
of  EAR  simulation  studies  to  include  (a)  design  of  experiments  and  (b)  an 
investigation  of  command  guidance  requirements. 

In  the  fourth-year  effort  presently  ending,  the  goal  of  the  first 
task  was  to  help  the  Army  in  the  selection  of  an  adaptive  tracking  filter 
which  could  track  a maneuvering  target.  The  goals  of  the  second  task  were 
to  investigate  quantization  effects  for  various  types  and  realizations  of 
digital  MTI  filter  and  to  recommend  a filter  which  represents  the  best  com- 
promise between  performance  and  hardware  requirements. 

Huntsville,  Alabama 
July  1975 
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ABSTRACT 

This  report  is  concerned  with  the  investiga  :ion  of  the  g-h-k 
filter  for  tracking  maneuvering  targets.  The  selection  of  the  filter 
coefficients  is  based  on  the  amounts  of  noise  and  maneuver,  and  other 
system  considerations  such  as  critical  damping  or  best  transient  response 
for  specified  smoothing.  Several  filter  initialization  schemes  were  tested. 
For  low  acceleration  maneuvers,  a considerable  amount  of  smoothing  can  be 
achieved  without  losing  track.  However,  in  order  to  track  severely  maneu- 
vering targets,  one  must  select  coefficients  which  give  a faster  transient 
response  at  the  expense  of  smoothing  capability.  Therefore,  it  is  logical 
to  use  an  adaptive  filter  with  a good  smoothing  capability  when  the  target 
is  not  maneuvering  and  a fast  response  during  a target  maneuver.  Clearly, 
the  main  problem  is  to  detect  the  maneuver  in  a reasonable  amount  of  time. 
This  can  be  done  using  a simplified  matched  filter  followed  by  a threshold 
detector.  The  proposed  adaptive  filter  was  evaluated  through  computer  simu- 
lation using  typical  trajectories.  The  performance  of  the  adaptive  filter 
is  limited  by  the  number  of  samples  required  by  the  detection  filter  and 
could  probably  be  improved  using  a more  complex  maneuver  detection  filter. 
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CHAPTER  1 
INTRODUCTION 

1.1  Prologue 

In  radar  tracking  systems,  the  raw  measurements  of  a target 
position  must  be  processed  to  provide  the  predicted  position  of  the  tar- 
get at  the  next  time  a radar  action  is  scheduled  for  this  target.  This 
estimate,  in  radar  coordinates,  is  then  used  to  steer  the  antenna  beam 
in  the  predicted  direction  at  the  next  tracking  instant.  In  order  to 
maintain  the  track,  the  prediction  accuracy  must  be  high  enough  that  the 
target  will  be  located  within  the  beam  and  within  the  range  bin  centered 
around  the  predicted  range.  Prediction  requires  not  only  a knowledge 
of  the  present  position,  which  can  be  measured,  but  also  an  estimate  of 
the  rate  (and  possibly  higher  derivatives)  of  the  target  motion,  which 
may  not  be  measured  by  the  radar.  In  view  of  the  fact  that  the  measure- 
ments are  obtained  at  discrete  intervals  of  time,  some  form  of  digital 
filter  suggests  itself  for  the  prediction  system. 

1.2  Tracking  Filters 

In  order  to  design  a tracking  filter,  it  is  necessary  to  have 
an  adequate  mathematical  model  for  the  motion  of  the  target.  In  practice, 
the  trajectory  will  not  be  a polynomial  for  all  time,  but  can  be  piece- 
wise  approximated  over  short  intervals  by  a polynomial  of  suitable  degree. 

Tracking  filters  [1]  may  be  non-recursive  (open-loop)  or  recursive  (closed- 
loop).  A well-known  filter  of  the  non-recursive  type  is  the  polynomial 
filter  12].  However,  the  most  widely  used  tracking  filters  are  of  the 
recursive  variety.  Prominent  among  these  are  the  Kalman  filter  [3,4], 
the  g-h  filter  [4-6]  and  the  g-h-k  filter  [7,8],  all  of  which  are  based 
on  polynomial  models.  The  Kalman  filter  minimizes  the  mean-squared  error 
if  the  model  is  exact  and  the  order  of  the  filter  is  one  more  than  the 
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degree  of  the  polynomial.  The  feedback  gains  of  a Kalman  filter  are 
different  at  each  step.  The  g-h  filter  is  a second-order  filter  and  is 
a simple  observer  [9]  for  constant  velocity  targets,  whereas  the  g-h-k 
filter  is  a third-order  observer  for  constant  acceleration  targets.  While 
the  gains  of  a Kalman  filter  are  derived  on  the  basis  of  statistical  con- 
siderations, those  of  the  g-h  and  g-h-k  filters  are  chosen  via  practical 
considerations.  In  their  most  common  forms,  the  g-h  and  g-h-k  filters 
have  constant  feedback  gains. 

The  design  of  a tracking  filter  must  take  into  account  two 
sources  of  error,  modeling  and  noise.  The  modeling  or  the  dynamic  error 
is  caused  by  the  target  motion  not  exactly  fitting  the  assumed  model.  This 
is  especially  severe  when  the  target  is  maneuvering,  either  due  to  atmo  - 
spheric  turbulence  or  due  to  an  evasive  maneuver.  The  noise  is  due  to  the 
radar  measurement  errors,  and  can  ba  adequately  modeled  by  a white  noise 
sequence.  The  requirements  for  rapid  dynamic  response  and  smoothing  of 
the  noise  are  always  in  conflict  and  a suitable  compromise  is  necessary  in 
designing  the  filter.  The  response  to  a changing  input  improves  as  the 
feedback  gains  are  made  larger,  but  this  results  in  a decreasing  ability 
to  smooth  out  random  fluctuations.  Other  factors  to  be  considered  in  de- 
signing a tracking  filter  include  the  choice  of  the  coordinate  system  and 
the  tracking  interval.  Among  the  various  coordinate  systems  available  [10] 
are  (i)  tne  spherical  coordinates,  range,  azimuth  and  elevation,  (ii)  the 
rectangular  coordinates  and  (iii)  a mixed  coordinate  system.  The  tracking 
interval  may  be  constant  or  varying. 

1.3  Scope  and  Organization  of  the  Report 

The  main  objective  of  the  task  reported  here  is  to  study  the 
applicability  of  g-h-k  filters  to  track  maneuvering  and  non-maneuvering 
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targets.  The  study  consists  of  two  parts:  (i)  the  selection  of  the  coeffi- 

cients g,  h and  k based  on  analytical  considerations  and  (ii)  the  develop- 
ment of  a procedure  to  dynamically  adjust  the  coefficients  to  adapt  to 
changes  in  t\e  trajectory. 

Throughout  this  report,  it  is  assumed  that  the  tracking  interval 
remains  constant,  «-hat  the  prediction  is  performed  in  the  range,  azimuth 
and  elevation  coordinates  and  that  the  three  filters  are  independent  of  each 
other.  It  will  also  be  assumed  that  only  position  measurements  and  no  ve- 
locity measurements  are  available. 

The  g-h-k  tracking  filter  is  briefly  reviewed  in  Chapter  2.  Input- 
output  relations  for  the  filter,  in  both  the  time-domain  and  the  frequency  «- 
domain  are  presented. 

The  selection  of  the  g-h-k  filter  coefficients  is  discussed  in 
Chapter  3.  Equations  are  derived  to  express  two  of  the  coefficients,  h and  k, 

in  terms  of  the  third  coefficient  g assuming  one  of  the  following  criteria: 

(i)  critical  damping,  (ii)  best  transient  response  for  specified  smoothing 
(optimal),  and  (iii)  use  of  steady-state  Kalman  gains.  The  choice  of  g should 
reflect  the  compromise  between  speed  of  response  and  smoothing  capability. 

The  performance  of  fixed  coefficient  g-h-k  filters  is  studied  in 

Chapter  4 via  simulation.  The  filters  are  applied  to  maneuvering  and  non- 
maneuvering  trajectories.  The  peak  and  RMS  errors  are  computed  for  dif- 
ferent maneuvers,  different  noise  strengths  and  for  various  values  of  the 
filter  parameter  g.  Both  optimal  filters  and  critically  damped  filters  are 
considered. 

A scheme  for  the  detection  of  a target  maneuver  is  proposed  in 
Chapter  5.  It  is  based  on  detecting  a bias  in  the  filter  residual  sequence 
and  is  an  approximation  to  the  decision  theoretical  method.  Also  presented 
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in  Chapter  5 is  a procedure  to  adjust  the  filter  coefficients  dynamically 
when  a maneuver  is  detected.  The  method  is  simulated  and  validated  using 
typical  trajectories. 

Chapter  6 contains  the  summary  and  conclusions  of  the  report 
along  with  recommendations  for  future  work. 


CHAPTER  2 
THE  G-H-K  FILTER 


Introduction 


The  g-h-k  filter,  also  known  as  the  a-$-y  filter,  is  a poly- 
nomial predictor-corrector  filter  of  the  third  order  and  is  of  the  linear 
recursive  type.  It  is  simply  an  observer  designed  to  reconstruct  the 
position,  velocity  and  acceleration  of  a constant  acceleration  target  based 
on  position  measurements.  Therefore,  it  can  track  a constant  acceleration 
trajectory  with  no  systematic  errors  in  the  steady  state.  In  addition  to 
the  predicted  estimate,  the  g-h-k  filter  also  provides  a smoothed  estimate 
of  the  present  position  which  may  be  useful  for  guidance  and  weapon  control 
functions.  In  its  standard  form,  the  g-h-k  filter  uses  only  the  position 
measurements,  and  must  be  modified  to  include  any  available  Doppler  measure- 


ments. 


2.2  The  Filter  Equations  [10] 

Let  x be,  in  general,  the  coordinate  in  which  the  filter  operates, 
so  that  x represents  either  range,  azimuth  or  elevation.  Assume  that  the 

smoothed  estimates  of  the  position,  velocity  and  acceleration  at  time  t^ 

? * 

are  known  and  denoted  as  2(n/n),  x(n/n)  and  x(n/n).  Then,  the  predicted 


state  at  time  tR+^  is  computed  as, 

x(n+l/n)  = x(n/n)  + T x(n/n)  + -7p  x(n/n) 

n z 

A A * 

x(n+l/n)  = x(n/n)  + T x(n/n) 

n 


(2.1) 


(2.2) 


x(n+l/n)  = x(n/n) 


(2.3) 


where 


T = t - t 
n n+1  n 


(2.4) 


After  the  position  measurement  x (n+1)  has  been  received  at  time  t the 

in  n+1 

predicted  state  can  be  corrected  resulting  in  the  smoothed  estimate  at  t 


£(n+l/n+l)  =*  £(n+l/n)  + g . . (x  (n+1)  - £(n+l/n)) 

n+x  ns 

x^n+l/n+1)  = x(n+l/n)  + (x  (n+1)  - x(n+l/n)) 

2kUl  “ 

x(n+l/n+l)  = x(n+l/n.)  + ---f'"  (x  (n+1 ) -x (n+1  ,/n) ) 

n+1 

where  gn+^»  ^n+i  an<*  ^n+l  are  t*ie  8ain  coefficients.  Denoting  the 
made  up  of  x,  x and  x by  X,  L.<L. , 


(2.5) 


(2.6) 


(2.7) 


vector 


X = x 


equations  (2.1)  to  (2.7)  can  be  written  as 


(2.8) 


X(n+l/n)  = $(n)  X(n/n) 


(2.9) 


X (n+1 /n+1)  = X(n+l/n)  + K(n+l)[x  (n+l)-H  X(n+l/n)] 
— — m — 


(2.10) 


where 


H = [ 1 0 0 ] 


and  the  state  transition  matrix  $(n)  is  given  by 


(2.11) 


(2.12) 


(2.13) 


Since  only  the  predicted  estimates  are  of  interest  here,  (2.9)  and  (2.10)  can 


be  combined  as 


where 


X(n+l/n)  = $(a) [I-K(n)H]  X(n/n-l)  + $(n)K(n)x  (n) 
— — m 

= ’i'(n)  X(n/n-l)  + <KnlK(n)x  (n) 

— m 


4\n)  = <^(r:)[I-K(u)H] 


(2.14) 

(2.15) 


and 


I is  the  third  order  identity  matrix. 


In  the  analysis  in  this  report,  it  will  be  assumed  that  the 


tracking  interval  T , and  the  gain  coefficients  g , h and  k are  constant, 
n n n n 

This  will  eliminate  the  dependency  of  ¥,  $ and  K on  the  time  variable  n, 
and  (2.15)  becomes, 


where 


2.3 


X(n+l/n)  = ¥ £(n/n-l)  + T x (n) 
— — m 


V = 


r = 


1 - g - h - k 

-(h+2k)/T 

-2k/T2 

g + h + k 
(h  + 2k) /T 
2k/T2 


T 

1 

0 


T2/2 

T 

1 


(2.17) 


(2.18) 


(2.19) 


Frequency  Domain  Representation 

The  system  represented  by  (2.17)  is  linear  and  shift  invariant 


with  x^(n)  as  the  input  sequence,  X(n+l/n)  as  the  vector  output  sequence. 
Therefore,  it  can  be  analyzed  in  the  frequency  domain  using  its  z-transfer 
function.  To  compute  the  transfer  function,  a z-transform  is  taken  on  both 
sides  of  (2.17)  which  yields, 

Z[X]  = z“1'FZ[X]  + rz[x  ] (2.20) 

~ — m 

where  the  operator  Z [ • ] denotes  the  z-transform.  Equation  (2.20)  can  be 
written  as, 

Z[X]  = (z  I-40"1  Zrz[x  ] . (2.21) 

Hi 

Using  (2.18)  and  (2.19),  the  above  equation  can  be  simplified  to  give  the 
transfer  functions  for  the  predicted  position,  velocity  and  acceleration. 


Each  transfer  function  is  of  the  type. 


H(z) 


3 2 

a^z  + a2z  + a^  z 
z 3 + b2z^  + b^z  + bg 


where  the  coefficients  are  given  in  Table  2.1. 


(2.22) 


TABLE  2.1  G-H-K  FILTER  TRANSFER  FUNCTION  COEFFICIENTS 


^^COEFFICIENTS 

a„ 

a. 

a. 

b„ 

b. 

b 

OUTPUT  TYPE^''~~--^^^ 

3 

2 

1 

2 

1 

0 

PREDICTED  POSITION 

5SH 

igggggggt 

g 

PREDICTED  VELOCITY 

h+2k 

T 

-2h-2k 

T 

ra 

g+h+k-3 

3-2g-h+k 

g-1 
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It  can  be  seen  that  the  characteristic  polynomial  (denominator  of  the 
transfer  function)  is  the  same  for  all  outputs  and  is  given  by 

c (z)  = z3- (3-g-h-k) z2+(3-2g-h+k) z~ (1-g)  . (2.23) 


Using  standard  z-domain  techniques  [11],  it  can  be  shown  that  the  system 
is  stable  if  the  following  conditions  are  satisfied: 

g > 
h > 
k > 

2g+h  < 

2g  > 
g (h+k)  > 2k  . 

It  can  also  be  shown  that  the  poles  of  the  transfer  function  move  towards 
the  unit  circle  as  the  coefficients  g,  h and  k are  decreased.  This  re- 
sults in  a reduced  bandwidth  which  provides  good  smoothing  at  the  same 
time  producing  a sluggish  dynamic  response.  On  the  other  hand,  an  increase 
in  the  filter  coefficients  improves  the  dynamic  response  at  the  cost  of 


0 

0 (2.24) 
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CHAPTER  3 

SELECTION  OF  THE  FILTER  COEFFICIENTS 

3.1  Introduction 

It  can  be  seen  from  the  filter  equations  presented  in  Chapter  2, 
that  the  g-h-k  filter  is  completely  defined  by  the  three  coefficients  g, 
h and  k.  As  was  discussed  earlier,  in  the  selection  of  the  coefficients, 
one  must  compromise  between  the  conflicting  requirements  of  good  noise 
smoothing  (narrow  bandwidth)  and  of  good  transient  capability  (wide  band- 
width). A very  practical  selection  scheme  would  be  to  express  two  of  the 
coefficients  in  terms  of  the  other  based  on  theoretical  considerations  and 
maintain  only  one  parameter  for  the  compromise.  The  choice  of  a value  for 
this  parameter  must  depend  on  the  system  application  where  t'-.e  relative  im- 
portance of  smoothing  and  dynamic  response  can  be  ascertained,  In  this 
chapter,  expressions  for  h and  k,  in  terms  of  g,  are  derived  based  on  three 
different  criteria.  The  merits  and  demerits  of  each  are  briefly  discussed. 
A fourth  coefficient  selection  scheme  is  described  which  should  be  useful 
during  heavy  transients. 

3 . 2 Critically  Damped  G-H-K  Filters 

Since  the  g-h-k  filcer  is  a third  order  observer  designed  to 
track  a constant  acceleration  trajectory,  the  steady-state  tracking  error 
for  such  a motion  is  zero  in  the  absence  of  noise.  However,  tracking 
errors  do  occur  during  transients  caused  by  improper  initialization  and  by 
changes  in  the  trajectory.  The  amount  of  time  required  for  the  transient 
errors  to  die  down,  is  a function  of  the  damping  of  the  filter.  An  under- 
damped filter  undergoes  periodic  oscillations  before  settling  down.  Under- 
damping is  caused  by  the  presence  of  a pair  of  complex  conjugate  poles  and 
can  be  avoided  by  forcing  all  the  poles  to  be  real.  Three  real  and  unequal 
poles  produce  no  oscillations  but  the  response  will  be  slow,  resulting  in 
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overdamping.  Critical  damping  provides  a fast  response  while  still  pre- 
venting oscillations. 

Critical  damping  of  a third  order  system  can  be  achieved  by 
forcing  all  the  poles  to  be  positive,  real  and  equal.  This  can  be  accom- 
plished by  selecting  the  coefficients  g,  h and  k such  that  the  character- 
istic equation  given  by  (2.23)  has  three  real  and  equal  roots.  Assuming 
the  roots  are  given  by  z = 8,  critical  damping  is  obtained  if  the  following 
equation  is  satisfied: 

(z— 3) 3 = z3- (3-g-h-k) z2+(3-2g-h+k) z- (1-g) . (3.1) 

Equating  the  coefficients  of  like  powers  of  z on  the  two  sides  of  (3.1) 
yields, 


3B  = 3-g-h-k 

(3.2) 

3B2  = 3-2g~h+k 

(3.3) 

and 

B3  - 1-g. 

(3.4) 

Solving  (3.2)— (3.4) 

for  g,  h and  k yields, 

g - i - e3 

(3.5) 

h = 1.5(1  - B2)(l-B) 

(3.6) 

k = 0.5  (1-B)3 

(3.7) 

Given  the  value  of  g,  (3.5)  can  be  solved  for  B which  is  then  used  to 
determine  the  values  of  h and  k from  (3.6)  and  (3.7).  Alternatively, 

B can  be  retained  as  the  compromising  parameter.  Stability  requires  that 


0 < & < 1 . (3.8) 

When  B = 0,  no  smoothing  is  applied  to  the  data  and  the  amount  of  smooth- 
ing increases  as  B increases.  When  3=1,  all  the  filter  coefficients 
are  zero  and  the  filter  output  becomes  independent  of  the  input.  Fig- 
ure 3.1  shows  the  variation  of  h and  k as  g varies  from  0 to  1. 
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The  stringent  condition  that  there  be  no  oscillations  could 
result  in  an  excessive  transient  period,  especially  following  a severe 
maneuver.  In  such  cases,  it  may  be  advantageous  to  allow  oome  oscilla- 
tions in  order  to  improve  the  dynamic  response. 

3.3  The  Optimal  Filter 

Here  again,  only  one  parameter  will  be  left  free  for  the  com- 
promise. But  for  every  value  of  this  parameter,  the  other  two  will  be 
chosen  to  give  the  maximum  noise  smoothing  for  a given  transient  capa- 
bility. It  is  in  this  sense  that  the  filter  is  optimal  and  the  synthesis 
procedure  is  similar  to  the  one  used  in  [6]  for  g-h  filters.  Before 
attempting  a compromise,  it  is  necessary  to  define  suitable  measures  of 
noise  reduction  and  transient  performance. 

For  noise  smoothing,  the  performance  measure  will  be  the  vari- 
ance reduction  ratio,  £,&. , the  steady-state  variance  in  position  output 
when  the  input  is  a unit  variance  white  noise  sequence.  If  { g (k) } Is  the 
impulse  response  sequence  for  the  predicted  position,  then  the  output  for 

an  input  sequence  of  { f (k) } is 

k 

c(k)  = l g(j)  f(k-j)  • (3.9) 

3=0 

The  autocorrelation  of  the  output  is, 

K(n)  .=  E[c(k)c(k+n)] 
k k+n 

= l l g(i)g(j)E[f (k-j)f (k+n-i)]  (3.10) 

3=0  i=0 

where  E[ • ] denotes  the  statistical  expectation  operator.  If  the  input 

is  a zero  mean  independent  sequence,  (3.10)  becomes,  in  steady  state, 

00 

K(n)  = l g(j)g(j+n)  . 

3=0 


(3.11) 
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Thus,  the  variance  of  the  output  sequence  is  given  by, 

00 

K(0)  = I g2(j)  . (3.12) 

j=0 

For  transient  performance,  the  measure  is  based  on  the  capa- 
bility to  follow  a constant  acceleration  input.  The  input  considered  is, 

x(k)  = (k2  + k)/2  . (3.13) 

If  the  predicted  position  output  sequence  is  (r(k)},  the  transient  per- 
formance measure  is, 

00 

D = l [x(k)  - r(k-l)]2  . (3.14) 

k=0 

This  simple  input  is  chosen  to  facilitate  analysis,  however,  it  retains 
the  important  feature  of  constant  acceleration. 

The  optimum  filter  is  synthesized  by  selecting  the  impulse  re- 
sponse sequence  { g (k) } so  as  to  minimize  K(0)  for  a given  D,  or  vice  versa. 
The  problem  can  be  solved  in  the  framework  of  constrained  optimization  in 
which  one  wishes  to  minimize 


J = K(0)  + XD 

where  X is  the  Lagrange  multiplier. 

The  z-transform  of  the  input  sequence  (3.13)  is, 

2 


(3.15) 


X(z)  = 


(z-1)' 


(3.16) 


Therefore,  if  G(z)  is  the  transfer  function  (z-transform  of  the  sequence 
{ g (k) } ) , the  z-transform  of  the  output  sequence  is  giver,  by 


R(z)  = X(z)  G(z) 

2 

— G(z) 


(3.17) 


Rearranging  (3.17)  yields, 

G(z)  = 


(z-1)' 


3 2 

z -3z  +3z-l 


R(z) 


(3.13) 
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An  Inverse  z-transfonn  of  (3.18)  gives, 

g(k)  = r(k+l)  -3r(k)  +3r(k-l)-r(k-2)  . (3.19) 

Using  (3.19),  (3.12)  and  (3.14),  (3.15)  becomes, 


J = l [ { r (k+1) -3r (k)+3r (k-1 ) -r (k-2 ) } 1 

k=0 


(3.20) 


+ X{x(k)-r(k)}  ] . 

In  order  to  find  the  first  variation  of  J,  the  optimal  r(k)  is  perturbed  by 

r (k)  -*•  r(k)  + e h(k) 

where  e is  a constant  and  h(k)  is  an  arbitrary  variation.  The  optimal 
r(k)  is  found  by  equating  the  first  variation  of  J to  zero,  L.Z., 


3J[r(n)  j £ h(n)] 
3e 


■ 0 


(3.21) 


0.22) 


for  k < 0 


(3.23) 


Substituting  (3.20)  in  (3.21)  yields, 

00 

l [ { r (k+1) -3r (k)+3r (k-1) -r (k-2 ) } { h (k+1) -3h (k)+3h (k-1) -h (k-2) } 
k=0 

+ X{x(k)-r(k-l)}(-h(k)}]=  0 . 

The  realizability  of  the  filter  requires  that, 

h(k)  = o' 
and  r(k)  - 0 

Using  (3.23)  and  rearranging  the  summation  indices,  (3.22)  becomes, 

00 

l h(k) [r (k+2)-6r(k+l)+15r(k)-20r (k-l)+15r(k-2) 
k=0 

-6r(k-3)+r(k-4)-X{x(k)-r (k-1)}]  = 0 . (3.24) 

Since  the  variation  h(k)  is  arbitrary,  the  quantity  inside  the  brackets 
must  be  identically  zero,  -t.e., 

r (k+2 ) -6r (k+1) vl5r (k)-20r (k-l)+15r (k-2) -6r (k-3) 

+r(k-4)+Xr (k-1)  - Xx(k)  = 0 (3.25) 


for  all  k > 0. 
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Because  the  left  hand  side  of  (3.25)  is  zero  for  all  k,  its  z-transform  can 
have  poles  only  on  or  outside  the  unit  circle.  Thus, 

F (z)  = (z2-6z+15-20z“1+15z_2-6z"3+1+Xz"1)R(z)-Xx(z)  (3.26) 

has  all  its  poles  on  or  outside  the  unit  circle.  Rewriting  (3.26),  yields 

F(z)  = + ..  R(z)  - X X(z).  (3.27) 

z 

Since  X(z)  given  by  (3.16)  has  no  poles  inside  the  unit  circle,  the  real- 
izable poles  of  R(z)  must  match  the  roots  of 

(z-1)6  + X z3  =0  . (3.28) 

It  can  be  easily  seen  that  if  z = a is  a root  of  (3.28),  then  z = 1/a  is 
also  a root.  Let  a,  b and  c be  the  three  realizable  poles  of  G(z), 

-i.C.. , a,  b and  c are  the  roots  of  the  characteristic  equation 


z3- (3-g-h-k) z2+(3-2g-h+k) z- (1-g)  = 0 . (3.29) 

Since  the  realizable  poles  of  R(z)  and  G(z)  are  th;  same,  it  follows  that, 
(z-a) (z-b) (z-c) (z-l/a) (z-l/b) (z-l/c)  = (z~l)6-Xz3  (3.30) 

where 

a+b+c  = 3-g-h-k  (3.31) 

a b c = 1 - g (3.32) 

ab  + be  + ac  = 3 - 2g  - h + k , (3.33) 


Equating  coefficients  of  like  powers  of  z on  both  sides  of  (3,30)  and 
using  (3.31)  - (3.33),  yields 

g2  = h(2-g)  - gk  (3.34) 

and  h2  = k2  +4  gk.  (3.35) 

The  above  equations  can  be  solved  for  h and  k in  terms  of  g,  giving 
h - (2g3-4g2)  + >4g6-64g5+64g4 

8 (1-g) 


(3.36) 
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and  k = h<2"S)  - &2..  . (3.37) 

g 

This  method  of  selecting  the  coefficients  was  based  on  minimizing  K(0) 
for  a given  D.  Equations  (3.36)  and  (3.37)  are  the  expressions  for  h and 
k in  terms  of  g.  The  optimal  g depends  on  X,  which  is  selected  in  such  a 
way  that  the  value  of  D equals  the  specified  value.  Figure  3.1  shows  the 
variation  of  h and  k as  g varies  from  0 to  1.  Once  again,  high  values  of 
g result  in  poor  smoothing  and  the  amount  of  noise  reduction  increases  as 
g decreases. 


This  set  of  coefficients  can  be  shown  to  result  in  a slightly 
underdamped  system.  The  optimality  of  the  system  depends  on  the  validity 
of  the  trajectory  model.  Therefore,  the  filter  is  very  well  suited  in  situa- 
tions where  the  motion  is  sufficiently  steady.  However,  if  the  target  accel- 
eration is  changing  rapidly,  its  transient  performance  may  not  be  adequate. 
3.4  Use  of  Steady  State  Kalman  Gains 

In  the  case  of  the  optimum  filter  described  in  the  previous  sec- 
tion, the  squared  error  in  the  predicted  position  was  minimized  assuming 
a deterministic  model.  The  Kalman  filter  [12]  is  based  on  a statistical 
model  and  minimizes  the  mean  squared  errors  in  the  three  outputs. 

The  message  r del  under  consideration  is  given  by 

-1)  * X(n)  + W(n)  (3.38) 

where  the  vector  X is  made  up  of  the  position,  velocity  and  acceleration 
as  defined  in  (2.8)  and  {w(n)l  is  an  independent  vector  noise  sequence  in- 
cluded to  account  for  modeling  inaccuracies.  The  state  transition  matrix 
$ is 


4>  = 


1 

0 

0 


T T2/2 

1 T 

0 1 


(3.39) 
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The  noise  sequence  (W(n)}  is  assumed  to  have  zero  mean  and  covariance 
matrix  Q,  i. £. , 


E[W(n) ] = 0 

E[W(k)WT(n)]  = Q 
= 0 


if  n = k 
otherwise  , 


(3.40) 


where 

^11 

Q12 

Q13 

Q = 

Q12 

Q22 

Q23 

Jl3 

Q23 

Q33_. 

The  measurement 

sequence 

x(n) 

m 

is  given 

by 

x (n)  = H X(n)  + v(n) 
m — 


(3.41) 

(3.42) 


where  H = [1  0 0] 

and  v(n)  is  a zero  mean  independent  noise  sequence  with  variance  R. 
Since  the  Kalman  filter  algorithms  have  been  widely  published,  only  the 
final  equations  will  be  repeated  here. 

' The  algorithm  for  the  predicted  estimate  at  the  (n+1)  th  stage 

based  on  measurements  up  to  and  including  the  n th  stage  is 

X(n+l/n)  = $>[  I-K(n)R]X(n/n-l)+4>K(n)xm(n)  . (3.43) 

.he  gain  vector  K(n)  is  computed  as 


m 


K(n)  = P(n)  HT[H  P(n)HT  + R]”1 
where  P(n)  is  the  predicted  error  covariance, 

P(n)  = E[(X(n)  - X(n/n-l) ){X(n)-X(n/n-l) }T  ] . 


(3.44) 


(3.45) 


The  error  covariance  matrix  is  computed  from  the  recursive  relation, 

P(n+1)  = $[I-K(n)H]P(n)  <J>T  + Q.  (3.46) 

A casual  examination  of  equations  (2.14)  and  (3,46)  reveals  that  the 
Kalman  filter  is  identical  to  the  g-h-k  filter  except  for  the  way  in 
which  the  gains  are  selected.  The  Kalman  filter  gain  K(n)  is  time  vary- 
ing as  can  be  seen  from  (3.44)  and  (3.46).  However,  it  can  be  shown  that, 
after  a large  number  of  steps,  the  error  covariance  matrix  and  hence  the 
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gain  vector  reach  steady  states  after  which  their  values  remain  constant. 
Therefore,  this  set  of  constant  gains  can  be  used  as  the  feedback  coeffi- 
cients of  the  g-h-k  filter. 

Let  and  denote  respectively,  the  elements  of  the  co- 
variance  matrix  and  the  gain  vector,  -i.e., 


!Pu(n) 

P12(n) 

Pl3(n) 

i 

P(n)  = 

;p12w 

P22<"> 

p23(n) 

(3  .47) 

.P13<n> 

P23(n) 

P33(n)- 

[h<*\ 

and 

K(n)  = 

K2(n)j  . 

(3.48) 

KjGO! 

Equation  (3.44)  can  now  be  expanded  to  give 


P11 (n) 

Vn)  ~ Pu(n)+  R 


4 -I  o 

K2(n)  = Pn(n)+  R 
P13(n) 

K3(n)  Pn(n)  + R 


(3.49) 


(3.50) 


(3.51) 


Using  (3.49)  - (3.51),  (3.46)  can  be  expanded  to  yield  [8], 

Pu(n+1)  = Kx(n)R  + 2 K2(n)  R T + K3  R T2  - P12(n)  K2(n)  T2 

- P13(n)  K2(n)  T3  + P22(n)  T2  + P?3(n)  T3 

- 1/4  P13(n)  K3(n)  T4  + 1/4  P3;J(n)  T4  + Qn  (3.52) 

P12(n+1)  + K2(n)  R + K3(n)  R T - p12(n)  K2(n)  T + P22(n)  T 
+ 3/2  P23(n)  T2  - 3/2  P12(n)  K3(n)  T2 

- 1/2  P13(n)  K3(n)  T3  + 1/2  P33(n)  T3  + Q12  (3.53) 

P13(n+1)  = K3(n)  P.  - ?12  (n)  K3(n)  T + P23(n)  T + 1/2  P33(n)  T2 

- 1/2  Pn(n)  K3(n)  T2  + Q13 


(3.54) 
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P22(n+1)  = -P12(n)  K2(n)  + P22(n)  + 2 P23(n)  T -2P12(n)  K3(n)  T 

-P13(n)  K3(n)  T2  + P^n)  T2  + Qn  (3.55) 

P23(n+1)  = -P12(n)  K3(n)  + P23(n)  - P13(n)  K3(n)  T + P33<n)T  + Q23  (3.56) 
P33(n+1)  = -P13(n)  K3(n)  + P33(n)  + Q33  . (3.57) 

In  order  to  find  the  steady  state  gains,  it  is  necessary  to  assign  values 
to  the  system  covariance  coefficients  Q^.  Assuming  that  the  target  is 
moving  with  constant  acceleration  except  for  an  uncertainty  in  the  acceler- 
ation, the  covariance  matrix  Q can  be  adequately  modelled  as. 


Q = 


0 0 
0 0 


0 

0 


0 


0 


(3.58) 


L m J 

2 

The  quantity  0^  is  the  maneuver  variance,  and  may  be  approximated  from  the 
maximum  expected  target  jerk.  Using  Equation  (3.58)  and  the  fact  that,  in 
steady  state, 


Equations  (3.49) 
gains  satisfy, 


and 


P^tn+l)  = P^  (n)  for  i,  j - 1,  2,  3 , 

- (3.57)  can  be  manipulated  to  show  that  the  steady  state 


2 T K3  = + K2  T + K3  Y 


2K1  K3 


(3.59) 

(3.60) 

(3.61) 


The  first  two  equations  can  be  used  to  express  two  coefficients  in  terms 
of  the  third,  which  is  then  chosen  using  the  last  equation.  Comparing 
(3.48)  with  the  g-h-k  filter  coefficients,  it  can  be  seen  that 


K, 


K„ 


h/T 

2k 


(3.62) 
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Use  of  (3.50)  in  (3.59)  and  (3.60)  yields, 

2h  = g(g  + h + k)  (3.63) 

h2  = 4 g k (3.64) 

which  can  be  solved  to  give, 

h = -2g  + 4-4  /l- g (3.65) 

k = h2/4g.  (3.66) 

Since  the  main  intent  of  this  analysis  is  to  express  two  coefficients 
in  terms  of  the  third,  (3.61)  will  not  be  discussed  further.  Given  a 
value  of  g,  (3.53)  and  (3.54)  can  be  used  to  compute  the  corresponding 
values  of  h and  k.  Figure  3.1  shows  the  variation  of  h and  k as  g varies 
from  0 to  1. 

As  can  be  seen,  the  steady  state  Kalman  gains  and  the  optimum 
gains  are  identical  except  at  large  values  of  g.  Since  the  steady  state 
gains  are  used,  it  should  be  expected  that  the  filter  has  a poor  transient 
performance.  Further,  the  optimality  of  the  Kalman  filter  depends  entirely 
on  the  validity  of  the  message  model.  Since  the  assumed  model  is  not 
adequate  when  there  is  a severe  maneuver,  the  filter  will  not  perform 
satisfactorily.  The  coefficient  selection  is  useful  mainly  during  steady 
periods,  where  high  noise  smoothing,  or  a small  value  of  g,  is  desired. 
However,  since  the  gains  are  the  same  as  the  optimum  gains  for  small  values 
of  g,  the  optimum  gains  could  be  used  in  all  situations. 

3 , 5 Coefficients  to  Improve  Transient  Performance 

All  the  three  selection  schemes  discussed  so  far  were  based  on 
optimizing  the  steady-state  performance.  However,  in  cases  where  heavy 
transients  are  expected,  it  may  be  desirable  to  design  the  filter  to  pro- 
vide a good  transient  response. 
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The  method  suggested  here  assumes  the  same  message  model, 

(3.38)  - (3.42),  as  the  Kalman  filter.  Using  this  model  and  the  equation 
(2.17)  of  the  g-h-k  filter,  it  can  be  shown,  by  elementary  linear  system 
covariance  analysis  methods,  that  the  error  covariance  matrix  P follows  the 
difference  equation, 

P(nfl)  = ¥P(n)  'i'1  + TR  rT  + Q • (3.67) 


The  covariance  matrix  P is  symmetric  and  hence  has  only  six  unknowns.  Let 
denote  the  vector  made  up  of  these  six  elements  of  P,  -L.Z., 

Pu(n)  - 


S(n)  - 


P12<n> 

P13(n> 

P22(n) 

P2.(n) 

P33(n) 


(3.68) 


where 
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Y = R 


and 


g 

gh/T 
2gk/T2 

hV 

2kh/T3 
4k2/T4 

Q11 
Q12 

Q13 

Q22 

q23 

Q33 

The  method  suggested  here  minimizes  the  cost  function 
N 

J - I Pu(n) 


o = 


n=0 


= l A S(n) 
n=0 


0 


0 


0 ] 


(3.72) 


(3.73) 


(3.74) 

(3.75) 

(3.76) 


where  A = fl  0 0 

In  (3.76),  N is  the  number  of  steps  over  which  the  transient  performance  is 
optimized.  A small  N results  in  improved  transient  performance,  whereas 
the  filter  reduces  to  the  one  described  in  Section  3.4  as  N tends  to  infinity. 
Assuming  the  initial  condition  S(0)  is  given  as  S!o,  the  minimization  can  be 
approached  from  an  optimal  control  point  of  view.  The  Hamiltonian  is  given 
by  H = A S(n)  - Ar(n+1)  [B  C S(n)  j Ay  + a)  (3.77) 

where  A is  the  Lagrange  multiplier  vector.  Tne  optimal  solution  satisfies, 


S(n+1)  = B C S(n)  + By  + a,  S(0)  = So 
T T 

A(n)  = A - C B A (n+1)  , A(N+1)  = 0 


(3.78) 


(3.79) 
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The  optimization  requires  the  solution  of  a two  point  boundary  value  prob- 
lem. Analytical  expressions  for  g,  h and  k cannot  be  obtained  since  the 
matrix  C and  the  vector  a contain  nonlinear  functions  of  g,  h and  k.  A 
simple  steepest  ascent  algorithm  similar  to  the  one  used  in  [13]  can  be 
utilized  to  solve  the  problem. 

3.6  Conclusions 

Three  methods  were  presented  where  two  of  the  filter  coefficients 
are  expressed  in  terms  of  the  third  in  such  a way  that  meaningful  criteria 
are  satisfied.  The  selection  of  the  third  coefficient  is  entirely  dependent 
on  the  practical  situation.  All  three  schemes  were  derived  using  a constant 
acceleration  trajectory.  Assuming  that  the  coefficients  remained  constant 
for  all  time,  the  methods  attempted  to  optimize  the  performance  for  all  time. 
However,  in  some  cases  where  heavy  transients  are  present,  it  may  be  de- 
sirable to  optimize  the  performance  over  a short  period  of  time.  A coeffi- 
cient selection  scheme  was  presented  where  this  is  accomplished  using  opti- 
mal control  principles. 
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CHAPTER  4 

PERFORMANCE  OF  THE  FIXED  COEFFICIENT  G-H-K  FILTER 

4.1  Introduction 

As  previously  discussed,  the  choice  of  a value  for  g is 
based  on  a compromise  between  speed  of  response  and  noise  smoothing  capa- 
bility. Another  factor  of  importance  is  the  selection  of  the  initial  state 
which  is  required  to  start  the  recursive  filter.  In  order  to  evaluate  the 
performance  of  the  fixed  coefficient  g-h-k  filter,  it  was  simulated  on  a digi- 
tal computer  and  applied  to  simulated  measurements.  The  purpose  of  this  ex- 
periment is  twofold:  (i)  to  select  a filter  initialization  scheme  which 

results  in  the  best  transient  performance  and  (ii)  to  compare  the  perform- 
ance of  the  filter,  using  several  values  of  g and  different  coefficient  se- 
lection schemes,  for  tracking  trajectories  with  varying  amounts  of  noise 
and  maneuver.  It  is  assumed  throughout  that  the  tracking  interval  is  con- 
stant and  equal  to  0.1  sec.  All  the  computer  programs  developed  for  this 
study  are  listed  in  Appendix  A. 

4.2  Measurement  and  Filter  Simulations 

4.2.1  Subroutine  TRAJ 

This  program  is  used  to  generate  samples  of  range  and  azimuth  of 
a maneuvering  or  nonmaneuvering  target.  It  is  assumed  that  the  target 
velocity,  the  tracking  interval  and  the  target  altitude  are  constant.  The 
simulated  trajectory  in  the  x-y  plane  consists  of  three  segments.  • a 
straight  line  path  from  time  t = 0 to  t = T^  seconds,  (ii)  a maneuver  to 
the  right  from  t = T^  to  t = T^  at  a constant  acceleration  and  (iii)  a 
straight  line  path  from  t = T^  to  t = T^.  Three  typical  trajectories 
which  are  used  in  the  analysis  are  shown  in  Figure  4.1.  Note  that  a non- 
maneuvering trajectory  can  be  obtained  by  setting  = T3.  In  addi- 

tion to  the  tracking  interval  DELT  which  is  transmitted  as  an  argument, 
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the  following  inputs  are  read  in: 

XO,  YO  the  initial  position  of  the  target 

Z the  constant  altitude 

VEL  the  constant  velocity 

THETA  the  initial  angle  measured  from  the  x-axis 

Tl,  T2,  T3  the  times  corresponding  to  T^,  T^  and  T^ 

AC  the  constant  radial  acceleration  during  the 

maneuver. 

The  outputs  of  the  program  are  the  number  of  samples  N,  and  two  arrays, 

R and  A,  consisting  of  N samples  each  of  the  range  and  azimuth.  All  the 
outputs  are  returned  through  arguments. 

The  equations  used  are: 
x(n)  = XO  + n.T.V  cos  6 

y (n)  = YO  + n.T.V  sin  0 
r = V^/a 


for  n T < Tl 


(4.1) 


a(n)  = V(n. t-Tl)/r 
Tl 

x(n)  = x(— )+r[cos(0-ir/2)-sin(a(n)+TT-0)  ] 
y(n)  = y(~^)+r[sin(0-ir/2)-cos(a(n)+/ir-0)  ] 


x(n)  = x(y)+(n.T-T2)V  cos(0-ct(y)) 


for  Tl<nT<T2 


/T2 


T2, 


y (n)  - y(-^)+(n.T-T2)V  sin(0-a(^)) 


for  T2<nT<T3 


(4.2) 


(4.3) 


where 


and 


T = DELT 
0 = THETA 
V = VEL 
a = AC 

r = radius  during  maneuver. 


The  range  and  azimuth  samples  are  calculated  as, 


R(n)  = / x2(n)  + y2 (n)  + Z2 
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(4.4) 

A(n)  = tan  1 (y(n)/x(n)).  (4.5) 

4.2.2  Subroutine  NOISE 

This  program  computes  the  simulated  measurements  by  adding  noise 
samples  to  the  range  and  azimuth  samples  generated  by  TRAJ.  The  noise  sam- 
ples form  a zero  mean,  Gaussian  independent  sequence  with  standard  devia- 
tion RSTD  and  ASTD  for  range  and  azimuth  respectively.  In  order  to  achieve 
a Monte-Carlo  type  simulation,  this  program  is  written  so  that  it  provides 
a different  noise  sequence  every  time.  In  addition  to  RSTD  and  ASTD,  the 
inputs  include  the  arrays  R and  AZ  which  contain  samples  of  range  and  azi- 
muth and  N the  number  of  samples.  The  outputs  are  the  arrays  XMR  and  XMA 
consisting  of  N samples  of  the  simulated  measurements  in  range  and  azimuth. 
All  inputs  and  outputs  are  transmitted  as  arguments.  The  noise  samples  are 
generated  using  the  subroutine  RANDN  of  the  UNIVAC  1108  library. 

4.2.3  Subroutine  GHKFIL 

This  program  simulates  the  g-h-k  filter  in  one  coordinate  based 
on  the  equations  given  in  Chapter  2.  An  input  INIT  is  used  to  select  one 
of  three  different  initialization  schemes  described  in  the  next  section. 
Other  inputs  are  the  array  of  measurements  XM,  the  number  of  samples  N, 
the  tracking  interval  DELT,  and  the  filter  coefficients  G,  H and  K.  The 
outputs  include  the  arrays  XP,  VP  and  AP  containing  predicted  position, 
velocity  and  acceleration,  respectively,  and  an  ar-**ay  RES  containing  the 
residual  sequence 

A 

RES(n)  = x (n)  - H X(n/n-l).  (4.6) 

m ~ 

All  inputs  and  outputs  are  transmitted  as  arguments. 

4.3  Filter  Initialization 


As  can  be  seen  from  the  equations  in  Chapter  2,  the  implementa- 
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tion  of  the  recursive  filter  requires  that  the  initial  states  be  specified. 
It  is  important  to  choose  these  states  properly  since  it  is  undesirable  to 
have  large  transients.  Assuming  that  three  measurements  have  been  made, 
three  different  initialization  schemes  are  considered: 


x(3/3)  = 


x (3) 
m 


(i) 


(2) 


5(3/3)  - 


x (3)  - x (2) 
m ra 


x(3/3) 

5(3/3) 


= 0 


5(3/3)  = 


5(3/3)  - 


x (3) 
m 

x (3)  - x (2) 
m m 


x (3)  - 2x  (2)+x  (1) 
m mm 


(4.7) 


(4.8) 


(3) 


x(3/3)  = 


5x  (3)+2x(2)-x(l) 
m mm 


A x (3)-x (1) 

x (3/3)  = ^ 


(4.9) 


x(3/3)  = 0 

The  first  set  of  equations  is  derived  assuming  a straight  line  between  the 
second  and  third  samples.  In  the  second  scheme,  the  position  and  velocity 
are  based  on  a straight  line  between  the  second  and  third  samples  while 
the  acceleration  assumes  piecewise  straight  lines  between  adjacent  samples. 

The  third  alternative  is  derived  through  a least  squares  straight  line  fit 
between  the  three  measurements.  By  setting  INIT  = 1,2,  or  3,  while  calling 
subroutine  GHKFIL,  the  corresponding  initialization  procedure  can  be  utilized. 

In  order  to  compare  the  three  initialization  schemes,  a main  pro- 


gram TSTINI  was  written  and  exercised  on  a nonmaneuvering  trajectory.  Since 
the  simulation  results  depend  upon  the  particular  noise  sequence  employed, 
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three  different  noise  sequences  were  used.  The  basis  of  comparison  was 
the  absolute  maximum  value  of  the  prediction  error  in  the  first  fifteen 
tracking  instances,  both  in  range  and  angle.  The  results  of  this  test 
are  given  in  Table  4.1.  The  nonmaneuvering  trajectory  shown  in  Figure  4.1 
was  used  along  with  the  following  parameters: 

RSTD  = 6m,  ASTD  = 7.5  m.rad 
g = 0.5  with  optimal  selection  of  h and  k. 

From  the  table  it  is  apparent  that  initialization  scheme  2 has 
the  worst  transient  performance.  This  is  probably  because  three  samples 
are  not  sufficient  to  estimate  the  acceleration.  Schemes  1 and  3 yield 
comparable  results  with  scheme  3 being  slightly  better.  Therefore,  it  will 
be  used  for  initialization  in  the  rest  of  the  analysis. 

4 . 4 G-H-K  Filter  Performance 

In  an  attempt  to  evaluate  the  performance  of  the  filter  under  a 
variety  of  conditions,  a main  program  TSTGHK  was  written.  The  program  was 
exercised  for  tracking  three  different  trajectories,  with  two  different 
noise  levels.  Several  values  of  the  coefficient  g were  tested  and  the 
coefficients  h and  k were  selected  using  both  the  optimal  and  the  critical 
damping  methods.  The  bases  of  comparison  are  the  absolute  maximum  error 
and  the  RMS  error  in  both  range  and  angle.  Once  again,  the  experiments 
were  repealed  for  three  different  noise  sequences.  Typical  results  are  given 
in  Tables  4.2  through  4.7  where  c and  0Q  denote  RSTD  and  ASTD  respectively. 
Since  noise  sequence  No.  3 provides  a performance  which  is  in  between  those 
of  the  other  two  sequences,  all  comparisons  are  made  with  respect  to  that 
sequence.  For  ease  of  comparison,  all  the  results  using  noise  sequence 
No.  3 are  repeated  in  Table  4.8.  The  following  different  values  of  the 
parameters  are  used  in  the  experiments  whose  results  are  shown  in  the  tables. 
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g : 

h,k  : 

Trajectory  : 
Noise  Levels: 


0.3  0.5  and  0.7 

optimal  and  critical  damping  methods 

three  trajectories  shown  in  Figure  4.1 

2 

no  maneuver,  30  m/see  maneuver  and 
60  m/sec^  maneuver 


a = 4m  and  oa  = 5 m.rad 
R v 

<t  = 8m  and  cQ  = 10  m.rad 
R vj 


Examination  of  the  results  reveals  the  following  observations: 

(1)  Low  values  of  g provide  the  best  performance  when  the 
target  is  not  maneuvering.  But  a target  maneuver  has  a 
deteriorating  effect. 

(2)  Target  maneuver  has  no  effect  on  the  errors  at  high 
values  of  g.  However,  they  have  very  poor  noise 
smoothing. 

(3)  At  low  values  of  g,  the  optimal  filter  is  better  than 
the  critically  damped  filter,  whereas  the  opposite  is 
true  at  high  values  of  g.  This  is  because  the  optimal 
filter  becomes  highly  underdamped  at  high  values  of  g. 

(4)  For  the  particular  noise  and  maneuver  parameters  chosen, 

a moderate  value  of  g,  g = 0.5,  seems  to  provide 

the  best  compromise. 

From  the  results,  it  is  apparent  that  a single  value  of  g cannot  provide  a 
satisfactory  performance  for  all  types  of  trajectories.  A small  value  of 
g is  adequate  for  tracking  nonmaneuvering  targets  but  may  lose  track  in 
case  of  a maneuver.  On  the  other  hand,  a large  value  of  g makes  the  per- 
formance insensitive  to  target  maneuvers,  but  may  lose  track  because  of 
its  poor  noise  smoothing  capability. 
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4.5  Conclusions 

The  performance  of  fixed  coefficient  g-h-k  filters  for  tracking 
maneuvering  and  nonmaneuvering  targets  was  studied  via  simulation.  Three 
different  filter  initialization  schemes  were  evaluated  in  the  light  of 
good  transient  performance.  The  tracking  capability  of  the  filter  was 
studied  by  varying  the  filter  and  target  parameters.  The  main  conclusion 
is  that  a single  set  of  filter  coefficients  is  not  sufficient  to  provide 
an  adequate  performance  for  both  maneuvering  and  nonmaneuvering  targets. 
Thus,  an  adaptive  filter  suggests  itself,  whereby  the  coefficients  are 
altered  when  changes  in  the  trajectory  are  detected.  A method  of  making 
the  filter  adaptive  to  target  maneuvers  is  discussed  in  the  next  chapter. 
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CHAPTER  5 

ADAPTIVE  G-H-K  FILTER  FOR  TRACKING  MANEUVERING  TARGETS 
5.1  Introduction 

The  most  commonly  used  tracking  filters,  i.  £. , the  Kalman  fil- 
ter, the  g-h  filter,  the  g-h-k  filter,  etc.,  operate  satisfactorily  pro- 
vided the  dynamic  model  on  which  the  tracking  algorithm  is  based  is  a 
correct  representation  of  the  true  state  of  the  target.  Usually,  the 
model  is  based  on  the  assumption  that  the  target  follows  a constant  ve- 
locity or  a constant  acceleration  trajectory.  However,  such  a model  can- 
not adequately  describe  a maneuvering  target,  not  only  because  the  con- 
stant acceleration  assumption  is  invalid,  but  also  due  to  the  inherent 
uncertainty  in  the  time  at  which  the  maneuver  is  initiated.  Singer  [3] 
has  suggested  that  a maneuvering  target  is  well  modeled  by  a linear  accel- 
eration model  driven  by  random  noise  of  proper  variance,  and  hence,  a 
Kalman  filter  can  be  used.  This  filter  maintains  track  throughout  the 
maneuver,  but  suffers  significant  degradation  when  the  target  is  not  maneu- 
vering. This  is  analagous  to  the  conclusion  drawn  in  the  previous  chapter 
that  a single  set  of  g-h-k  filter  coefficients  cannot  provide  a satisfactory 
performance  for  both  maneuvering  and  nonmaneuvering  targets. 

One  approach  towards  resolving  this  problem  is  the  use  of  adaptive 
trackers,  wherein  the  filter  senses  deviations  of  the  trajectory  from  the 
model  and  modifies  the  algorithm  accordingly.  An  adaptive  ’>'alman  filter  is 
proposed  in  [14]  where  the  predicted  estimate  is  the  weighted  sum  of  two 
Kalman  predictors,  each  of  which  is  based  on  a different  model.  The  rela- 
tive weights  are  based  on  a likelihood  ratio  which  reflects  the  possible 
occurrence  of  a maneuver.  In  the  adaptive  estimator  described  in  [15]  and 
[16],  target  maneuvers  are  modeled  as  unknown  random  variables  which  are 
estimated  along  with  the  target  state.  The  unsupervised  tracking  scheme 
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given  in  [17]  is  a feedback  filter  where  the  gains  are  selected  to  force 
orthogonality  in  the  residual  sequence.  The  gains  themselves  are  computed 
using  a stochastic  approximation  algorithm.  A decision-directed  adaptive 
tracker  is  explained  in  [18],  where  a maneuver  is  detected  from  the  bias 
in  the  residual  sequence  of  a Kalman  filter.  The  filter  is  reinitialized 
and  a different  value  is  used  for  the  system  noise  covariance  when  a maneu- 
ver is  detected.  Linear  regression  filters  are  used  in  [19]  for  tracking 
maneuvering  aircraft  targets  in  cartesian  coordinates. 

In  this  chapter,  an  attempt  is  made  to  design  an  adaptive  g-h-k 
filter.  The  procedure  is  to  use  a fixed  coefficient  g-h-k  filter  with 
high  noise  smoothing  capability  when  the  target  is  not  maneuvering.  How- 
ever, when  a maneuver  is  detected,  the  coefficients  are  modified  so  as  to 
provide  a good  transient  performance.  Target  maneuvers  are  detected  by 
processing  the  residual  sequence  of  the  g-h-k  filter.  The  procedure  is 
simulated  and  exercised  on  maneuvering  and  nonmaneuvering  trajectories  to 
ascertain  its  validity.  All  the  programs  used  are  listed  in  Appendix  B. 

5.2  Maneuver  Detection 

5.2.1  Principle  of  the  Method 

Using  the  same  notation  as  in  Chapter  3,  a maneuvering  target  is 
modeled  by  the  system  equations 


X(n+1) 
x (n) 

m 


<J>  X(n)  + A(n) 
H X(n)  + v(n) 
1 T 

$ = 0 1 

0 0 

H = [ 1 0 


2 

T /2 
T 
1 

0 ] 


(5.1) 

(5.2) 


where 
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and  the  vector  forcing  function  A(n)  is  included  to  account  for  the  tar- 
get maneuver.  Assuming  that  the  maneuver  was  initiated  at  n = n where 

o 

n is  unknown , 
o 

A(n)  = 0 for  n < nQ  . (5.3) 

The  filter  equation  is  given  by, 

A A 

X(n+l/n)  = <5  X(n/n-l)  + $ K v(n)  (5.4) 

where  v(n)  is  the  residual  sequence, 

A, 

V(n)  = xm(n)  - H X(n/n-l)  , (5.5) 

and  the  gain  vector  K is  given  by 


K = 


g 

h/T 


(5.6) 


Consider  a nonmaneuvering  system  given  by, 

Y(n+1)  = *Y(n)  (5.7) 

>'ffi(n)  = H Y(n)  + v(n)  (5.8) 

with  the  same  initial  condition  as  that  of  (5.1).  The  corresponding  pre- 
diction equation  is, 

A 

Y(n+l/n)  = 4>  Y(n/n-l)  + $K  V]L(n)  (5.9) 


where  v (n)  = y (n)  - H Y(n/n-l).  (5.10) 

1 m — 

It  can  be  easily  seen  that, 

A A 

X(n/n-l)  = Y (n/n-1)  ) (5.11) 

) for  n < n 


and  v(n)  = V^(n)  ) 

If  the  Kalman  gains  were  used  for  the  feedback  gains  in  (5.8),  it  io  well 
known  that  the  residual  sequence  would  be  a zero  mean  white  noise  sequence. 
While  this  is  not  true  in  general  for  the  g-h-k  filter,  under  the  condition 
that  the  filter  is  in  steady  state,  it  can  be  assumed  that  the  sequence  v^(r.) 
is  a zero-mean  independent  sequence.  However,  this  property  is  dependent  on 
the  model  (5.7)  being  a true  representation  of  the  trajectory.  If  the  target 
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is  maneuvering,  the  residual  sequence  develops  a bias  which  reflects  the 
maneuver.  The  residual  sequence  v(n)  can  be  written  as 

(5.12) 

(5.13) 

Therefore,  detecting  a bias  reduces  to  the  followin'  hypothesis  test: 

H 


v(n)  = v^(n)  + b(n) 


where  b(n)  = 0 for  n < n . 

~ o 


no  maneuver  : V(n)  = V^(n) 


: maneuver 


v(n)  = v^(n)  + b(n) 


(5.14) 


where  v^(n)  is  an  independent  sequence.  That  is,  detecting  a maneuver  is 
equivalent  to  detecting  for  the  presence  of  a deterministic  signal  of  un- 
known amplitude  and  time  of  arrival  in  a background  of  zero-mean  white 
noise. 

5.2.2  Equation  for  the  Bias 

Comparison  of  (5.2)  and  (5.8)  gives 

(5.15) 

Combination  of  (5.5)  and  (5.15)  yields  , 

A 

v(n)  = y (n)  - H Y(n)  + H X(n)  - H X(n/n-l) 
m — — ~ 

A A A 

= ym(n)  - H Y(n/n-l)  + H(X(n)  - Y(n))  - H(X(n/n-l)-Y(n/n-l)) 


xm(n)  = ym(n)  - H Y(n)  4-  H X(n) 


- vx(n)  + H(X?(n)  - X3(n)) 


where 


X2(n)  = X(n)  - Y(n) 

A A 

X (n)  - X(n/n-l)  - Y(n/n-l)  . 

From  (5.1)  and  (5.7),  it  is  easy  to  show  that, 

X2(n+1)  - $X2(n)  + A(n) 

and  (5.4)  and  (5.9)  can  be  manipulated  to  give, 

X (n  i-1)  = *(j-KH)X,(n)  + H X (n) , 

’""w  4 

Defining, 


(5.16) 

(5.17) 

(5.18) 

(3.19) 

(5.20) 


X,  (nl  = X„  (n)  - X„(n) 

*+  jL 


(5.21) 
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equations  (5.19)  and  (5.20)  can  be  combined  to  yield, 

X^n+l)  = $(I-K  H)X^(n)  + A(n) . (5.21) 
Comparison  of  (5.12),  (5.16)  and  (5.21)  shows  that  the  bias  sequence  can 
be  expressed  as, 


b(n)  = H X^n+l)  . (5.23) 
Since  the  maneuvering  and  nonmaneuvering  models  assume  the  same  initial 
conditions, 

2^(0)  “ 0 . (5.24) 
Equations  (5.22),  (5.23)  and  (5.24)  provide  a model  for  the  bias  in  the 
residual  sequence. 

The  initial  effect  of  a target  maneuver  can  be  adequately 
modeled  by  a sudden  change  in  acceleration,  -i. £., 


A(n) 


0 

a 


for  n = n 


(5.25) 


= £ otherwise 

Using  (5.25),  it  can  be  shown  that  the  z-transform  of  the  bias  is  given 


fay»  _n 

B(z)  = g z °_lz±l_) (5.26) 

z3  -(3-g-h-k)z2  + (3-2g-h+k) z- (1-g) 

Assuming  that  the  filter  coefficients  are  selected  to  provide  critical 
damping  according  to  (3.5)  - (3.7),  the  bias  sequence  can  be  shown  to  be 
of  the  form, 

b (n)  = c T2(n-nQ)2  B(n''n°)  (5.27) 

where  c is  a constant.  A large  value  of  g results  in  a small  value  of  B 
and  hence  a small  bias.  This  is  in  agreement  with  the  earlier  result  that  a 
large  value  of  g produces  a fast  reaction  to  changes  in  the  trajectory. 

5.2.3  Bias  Detection 


As  stated  before,  the  maneuver  detection  reduces  to  a hypothesis 
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test  where  it  is  required  to  detect  the  presence  of  a deterministic  signal 

b(n)  and  its  time  of  occurrence  n in  a background  of  white  noi^e-  The 

o 

classical  decision  theoretical  method  is  to  pass  the  residual  sequence  v(n) 

through  a matched  filter  [32]  which  is  matched  to, 

m(n)  = (nT)2  gn  (5.28) 

and  subject  the  output  to  a threshold  detector.  Assuming  that  the  detection 

must  be  completed  within  p steps  of  the  occurrence  of  the  maneuver,  the 

impulse  response  of  the  matched  filter  is, 

y(n)  = m(p  - n)  0 < n £ p (5  29) 

= 0 otherwise  . 

It  is  desirable  to  implement  the  matched  filter  as  a recursive  filter  so 
as  to  avoid  excessive  computation  and  storage.  An  impulse  response  of  the 
type  given  in  (5.29),  requires  a filter  of  order  k.  In  order  to  simplify 
the  implementation,  the  impulse  response  y(n)  is  approximated  by  an  expo- 
nential, which  can  be  synthesized  as  a first  order  recursive  filter.  De- 
noting the  new  impulse  response  as, 

A(n)  = an  , (5.30) 

the  difference  equation  governing  the  matched  filter  is  given  by, 

p(n+l)  = a y(n)  +T  v(n+l)  (5.31) 

where  p(n)  is  the  output  and  the  input  is  the  residual  sequence  v(n).  The 
output  sequence  p(n)  is  passed  through  a threshold  detector.  If  the  out- 
put exceeds  the  threshold  at  n = £,  a maneuver  is  detected  and  the  time  of 
its  occurrence  is  estimated  as 

nQ  = l - p . (5.32) 

The  remaining  problems  are  the  selection  of  proper  values  for  p and  a. 

A large  p results  in  more  maneuver  data  in  the  matched  filter,  thereby 
making  it  easier  to  detect,  but  the  lag  between  the  occurance  of  the  maneuver 
and  its  detection  may  be  unacceptable.  If  p is  small,  detection  is  based  on 
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few  samples  containing  maneuver  information,  and  the  performance  could  be 
quite  poor.  The  value  of  a should  be  chosen  so  that  A(n)  approximates 
y(n)  and  is  therefore  a function  of  g. 

5.3  Adaptive  Filtering  Scheme 

The  maneuver  detection  method  explained  in  the  previous  section 
is  utilized  to  design  a g-h-k  filter  which  is  adaptive  to  target  maneuvers. 

It  is  assumed  that  the  target  is  initially  on  a straight  line,  path  and  that 
it  is  being  tracked  with  a g-h-k  filter  using  a small  value  for  g to  pro- 
vide good  noise  smoothing.  The  residual  sequence  v(n)  is  passed  through 
the  first  order  filter  given  by  (5.31)  for  maneuver  detection.  The  output 
y(n)  of  this  filter  is  compared  to  a preselected  threshold  6.  If  at  n = n^, 
y(n)  = y(n^)  exceeds  the  threshold  6,  a target  maneuver  is  indicated.  In 
such  a case,  the  filter  coefficients  are  increased  in  order  to  provide  a 
fast  response  to  trajectory  changes  and  the  filter  is  reinitialized.  In 

order  to  reinitialize  the  filter,  the  I latest  measurements  and  smoothed 

s 

values  are  stored.  When  a maneuver  is  detected,  the  filter  is  applied  using 
the  new  coefficients  and  starting  from  the  earliest  available  smoothed  values. 
Since  large  values  of  coefficients  also  result  in  poor  noise  smoothing,  the 
filter  coefficients  are  switched  back  to  the  original  values  after  1^  track- 
ing Intervals  following  the  maneuver  detection. 

Subroutine  GHKADA  implements  the  above  schem-  for  one  radar  coordi- 
nate. The  input  to  the  program  transferred  through  arguments  are  the  array 
of  measurements  XM,  and  the  threshold  6 for  detecting  the  maneuver,  denoted 
as  THR.  Other  inputs  transferred  through  a common  statement  are  the  number 
of  samples  N,  the  tracking  interval  DELT,  the  filter  coefficients  GN,  HN, 

KN  for  the  nonmaneuvering  part,  the  coefficients  GM,  HM,  KM  tor  the  maneuver- 
ing part,  the  detection  filter  weight  a (see  (5.31))  denoted  as  WEIGHT,  and 
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and  the  values  of  I and  1^  denoted  as  ISTORE  and  ITRANS,  respectively. 

The  outputs  are  the  arrays  XP,  XV  and  XA  containing  predicted  position, 
velocity  and  acceleration,  respectively,  and  are  transmitted  as  arguments. 

Main  program  TSTADA  is  similar  to  TSTGHK  developed  in  Chapter  4. 

It  is  used  to  find  the  maximum  and  the  RMS  values  of  the  error  in  the  pre- 
dicted position  for  various  trajectories  generated  using  subroutine  TRAJ. 

5.4  Adaptive  Filter  Performance 

To  evaluate  the  performance  of  the  adaptive  scheme,  the  sif’la- 
tion  program  was  exercised  on  typical  trajectories  described  in  Chapter  4. 
The  experiment  was  repeated  for  various  combinations  of  the  filter  coeffi- 
cients and  target  maneuvers.  A typical  set  of  results  is  presented  in 
Table  5.1.  During  the  nonmaneuvering  position  of  the  trajectory  the  op- 
timal coefficients  with  GN  = 0.3  are  used.  The  value  of  the  detection 
filter  weighting  coefficient  WEIGHT  was  chosen  as  0.85,  since  this  selection 
provides  a good  approximation  to  the  matched  filter  when  g = 0.3.  The 
threshold  for  maneuver  detection  in  range  was  set  at  6 = 4.0,  which  was 
selected  based  on  preliminary  experimental  results.  For  the  types  of  tra- 
jectory used,  the  target  maneuver  does  not  have  an  appreciable  effect  on 
the  angle  tracking  filter.  Therefore  only  the  results  of  the  range  track- 
ing filter  are  given.  The  table  gives  the  peak  range  error  with  GN  = 0.3 
using  the  adaptive  and  the  nonadaptive  filter.  The  adaptive  filter  uses 
GM  = 0.5  and  0.6  with  both  the  optimal  and  critical  damping  coefficient 
selection  schemes.  The  values  of  I and  I used  are  3 and  10,  respectively. 

e ^ 

The  results  of  Table  5.1  and  others  not  presented  permit  the 
following  inferences. 

(1)  The  adaptive  scheme  reduces  the  maximum  error  for  high 


acceleration  maneuvers. 


TABLE  5.1  PERFORMANCE  OF  ADAPTIVE  FILTER 
(PEAK  RANGE  ERROR  IN  METERS) 
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(2)  During  the  maneuver,  the  selection  of  the  filter 
coefficients  based  on  critical  damping  is  better  than 
than  the  optimal  selection.  This  is  to  be  expected 
since  the  optimal  scheme  is  slightly  underdamped  and 
therefore,  not  ideal  during  transients. 

(3)  The  selection  of  I is  not  critical.  In  fact,  the 

s 

simpler  choice  of  I = 1 gives  almost  an  identical 
performance  to  I = 3. 

(4)  The  maximum  error  occurs  before  the  maneuver  is  de- 
tected. Therefore,  the  adaptive  scheme  would  work 
better  if  the  maneuver  were  detected  earlier. 


50 


CHAPTER  6 

SUMMARY,  CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  EFFORTS 

6.1  Summary 

Chapter  2 is  a review  of  the  g-h-k  tracking  filter  and  presents 
all  the  relevant  equations.  The  selection  of  the  filter  coefficients  is 
discussed  in  Chapter  3.  Two  of  the  filter  coefficients  can  be  expressed  in 
terms  of  the  third  based  on  performance  criteria,  which  involve  noise  smooth- 
ing and  speed  of  response.  The  purpose  of  Chapter  4 is  twofold:  (i)  to 

select  a filter  initialization  scheme  for  best  transient  performance  and 
(ii)  to  select  the  filter  coefficients  for  various  amounts  of  noise  and 
maneuver.  Since  the  best  choice  of  the  filter  coefficients  for  nonmaneu- 
vering and  for  maneuvering  targets  is  quite  different,  it  seems  logical  to 
make  the  filter  adaptive,  that  is,  to  change  the  coefficients  as  soon  as 
a maneuver  is  detected.  Chapter  5 proposes  a technique  for  maneuver  detec- 
tion which  is  based  on  the  bias  in  the  residual  sequence.  The  performance 
of  the  adaptive  filter  is  evaluated  through  computer  simulation. 

6.2  Conclusions 

The  selection  of  the  g-h-k  filter  coefficients  is  very  important 
and  must  be  a compromise  between  the  conflicting  requirements  of  good 
noise  smoothing  and  of  good  transient  capability.  Assuming  either  critical 
damping  or  optimal  transient  response  for  a specified  amount  of  smoothing, 
one  can  express  h and  k in  terms  of  g,  where  g is  chosen  to  satisfy  the 
smoothing/ transient  compromise. 

The  performance  of  the  fixed  coefficient  g-h-k  filter  is  quite 
good  for  a nonmaneuvering  target  and  is  still  adequate  for  a wide  range 
of  maneuvers.  The  filter  initialization  has  a serious  effect  on  the  per- 


formance and  a scheme  is  proposed  based  on  a straight  line  fit  through  the 
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three  measurements.  At  low  values  of  g,  the  optimal  filter  is  better  than 
the  critically  damped  filter,  whereas  the  opposite  is  true  at  high  values 
of  g.  This  is  because  the  optimal  filter  becomes  very  underdamped  at  high 
values  of  g.  A large  value  of  g provides  the  best  performance  when  the 
target  is  maneuvering  and  should  be  used  in  conjunction  with  critical  damping. 

Trajectories  were  generated  where  a straight  line  path  is  fol- 
lowed by  a high  acceleration  maneuver.  It  was  shown  that  the  maneuver  in- 
troduces a bias  in  the  residual  sequence.  A technique  which  approximates 
matched  filtering  followed  by  threshold  detection  was  proposed  to  detect 
the  maneuver.  The  detection  requires  at  least  eight  samples  after  the 
maneuver  starts.  An  adpative  filtering  scheme  was  suggested  where  the 
coefficients  are  selected  for  good  smoothing  performance  prior  to  maneuver- 
ing and  for  good  transient  performance  following  a maneuver  detection.  It 
is  necessary  to  reinitialize  the  filter  when  the  coefficients  are  changed. 

The  use  of  the  adaptive  filter  reduces  the  tracking  error  during  maneuver- 
ing. The  maximum  error  occurs  before  detection  of  the  maneuver. 

6.3  Suggestions  for  Future  Efforts 

The  performance  of  the  adaptive  filter  is  limited  by  the  detec- 
tion capability.  It  could  be  greatly  improved  if  a maneuver  were  detected 
in  fewer  samples.  Therefore,  one  may  want  to  use  a better  approximation 
for  the  matched  filter.  The  fixed  threshold  was  determined  experimentally, 
and  a more  systematic  procedure  is  recommended. 

All  the  study  assumes  a constant  tracking  interval.  However,  if 
the  radar  ir  used  to  track  several  targets,  it  may  be  advantageous  to  in- 
crease the  sampling  rate  for  a short  interval  following  maneuver  detection. 
One  could  consider  other  reinitialization  techniques  than  the  one  described. 
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t. 


The  filter  was  made  adaptive  to  target  maneuvers.  In  addition, 
one  could  also  make  the  selection  of  the  g-h-k  coefficients  a function  of 
the  RMS  value  of  the  noise  to  be  estimated  from  the  residual  sequence. 


APPENDIX  A 


PROGRAMS  FOR  EVALUATION  OF  FIXED  COEFFICIENT  G-H-K  FILTERS 

.nis  appendix  contains  the  FORTRAN  listings  of  the  programs 
which  are  used  for  the  evaluation  of  g-li-k  filters.  The  programs  are: 
subroutines  TRAJ,  NOISE  and  GHKFIL,  and  main  programs  TSTINI  and  TSTGHK 
All  the  programs  are  described  in  Chapter  4. 
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SUBROUTINE  TRAJ(RiA,N,DELT) 

DIMENSION  RIIliMl) 

READ  1 »XO.YQ.Z,VEL, THETA 
PRINT  2,X0,Y0.Z.VEL,THETA 
READ  1 , T 1 ,T2,T3,AC 
PRINT  3 , T 1 iT2.T3.AC 

1 FORMAT ( ) 

2 FORMAT { 5X » * XO  YO  Z VEU  THETA  *,5F12.5) 

3 FORM AT ( 5X  » • Tl  T2  T3  AC  * • 4F 1 2. 5 ) 

T»0. 

N*U 

CT "COS (THETA ) 

ST«SIN(THETA) 

IFUC.GT.l.)  rad«vel»vel/ac 
IQ  IF(T.6T.T1 ) 60  TO  11 
N“N+  1 
VT-VEL*T 
X»XU*VT*CT 
Y«YO+VT*ST 

R(N)«SQRTCX**2*Y**2+Z**2) 

AIN) « AT AN2 ( Y # X ) 

T-T+DELT 
60  TO  10 

11  CONTINUE 

IF  * AC  «LT  • 1 • > RETURN 

X i »x 

Yl«Y 

12  IFIT.6T.T2I  60  TO  13 
alpha-velmt-ti  )/RAD 

N«N+  1 

AN6»ALPHA-THETA 
X«XI*RAD*IST*SINUN6)  I 
Y«Yl*RAD*(-CT«-C0S(AN6)» 
R(N)«SQRT(X**2*Y**2*Z**2> 
A(N)«ATAN2(Y,X) 

T-T+DELT 
60  TO  12 

13  CONTINUE 
X2«X 
Y2»Y 

phj-theta-alpha 
CP«COS(PHI ) 

SP»S I N ( PHI ) 

14  IFIT.6T.T3)  RETURN 
N »N  ♦ 1 

V T*VEL* ( T-T2  ) 

X»X2*VT*CP 

Y»Y2«-VT*SP 

RIN)»SQRTIX«*2+Y**2+Z**2) 

A l N ) « AT  AN2 < Y , X ) 

t»t+delt 

60  TO  14 
END 
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SUBROUTINE  N01SE'H,AZ.N,XMR,XMA,HSTD,ASTD) 
01  HENS  I ON  R( i 5 ,AZ( 1 ) ,XMR(  1 ) ,XMA(  1 ) ,XN(H00u) 
DATA  SAVE/1R76./ 

XN ( l )«SAVE 
N2"2*N 

CALL  KAN0N(XNiN2iU»L'iU0) 

00  10  I«l,N 

XMR  ( n»R(  I )*RSTD*XN(  I ) 

10  X M A I I ) ■ A Z I 1 1*AST0*XW( I ♦ N ) 

SAVE«ABS(  XN(N2  H*H3S76 

RETURN 

END 


SUBROUTINE  GHKF I L ( XM , N , XP  , VP  , Ap  , UELT  . G , H , K , R£S  , I N I T ) 
DIMENSION  XM( 1 ) ,XPIl) . VP(1) , API  1 » ,RESI  II 
REAL  K 

GO  TO ( 20 • 2 1 |22I tINIT 
2U  X S ■ X M ( 3 ) 

VS» ( XM I 3 ) *XM ( 2 1 I/DELT 
AS«0. 

GO  TO  11 

21  X Sb XM  < 3 ) 

VSbIXM(3)-XM(2) l/OELT 
ASb(XM(3)-2.*XM{2»*XM{1))/(0ELT**2) 

GO  TO  11 

22  XS»(5.*XM<3)*2.*XM(2)-XM(l)>/6. 

V5M  XMI 3)-XM(l) )/ ( 2.*DELT  ) 

A b « U • 

11  CONTINUE 

DO  12  l ■*! » N 

XPI I )*XS*DELT*VS  + AS*10ELT**2/2# ) 

VPI  I >«VS*DEIT*AS 
API  I I « AS 

RESI 1 JbXMI I l-XPI I 1 
XS»XP( I ) +G  *RE  S I I I 
VSbVPI  I ) +HES  I 1 MH/OELT 

12  AS*  AP  I l)*RESm*(2t*K/0ELT**2) 

RETURN 

END 


56 


MAIN  PROGRAM  TSTGHK 


0  I MtNS  I UN  RI2UU)  ,AZ(2UU)  ,XNK(20U) ,XMA<200) 
DIMENSION  HP  <200)  , RV12U0)  ,RA(2U0)  , RES  (200) 
DIMENSION  AP(2UU),Av(200),AA(20U),AES|2QO> 
REAL  K 

1 FORMAT!  ) 

2 FORMAT ( B X • • G H K » ,3F12.b) 

J FORMAT (bXt*  R S T u ASTD  *,  *F 1 2 » b ) 

A f ORMaT ( bX  i 1 I N I T * ,ib) 
b F0RMAT(/5X,*  RESULTS  FOR  INCISE 


* lb// 

* * F I 2 • b/ 

* » F l 2 » b/ 

* iF12.b/ 

* * F I 2 • b/ 
• »Fl2.b/ 

* i F I 2 • S > 


• 20X » * RANGE  mean  SWUAkED  ERROR  * 

* 2UX,»  RANGE  RMS  ERROR  « 

* 20X,*  RANGE  MAXIMUM  ERROR  * 

* 20X i * ANuLE  MEAN  SWUARED  ERROR  * 

• 20X,»  ANGLE  R M S ERROR  a 

• 20X,«  ANGLE  MAXIMUM  ERROR  * 

7 FORMAT ( 5X i • CELT  • * F 12.5) 

read  i.delt 

PRINT  7.DELT 

READ  liGtH.K 

PRINT  2.G.H.K 

READ  1 i RS  TO  i ASTD 

PRINT  3iKSTU,AST0 

READ  1 » 1 N I T 

PRINT  **,1NIT 

CALL  TRAJ»R,aZ|W.CELT) 

00  K-  I NOISE*  1,3 
S 5 '•)  R * G » 

SSC|A*0  « 

ERMAXR-U. 

ERMAXAaO, 

call  noise (k,aZ,n,Xmr ,xma,rsto,mSTD) 

call  GHK.FI L ( XmR,N ,RP,RV ,ra ,DELT  ,G,H|K , res,  ini  r ) 

CALL  GHnFIL(XMA,N,AP,AV,AA,DELT,G,H,K,AES,INIT) 

DO  11  1«U,N 

RE*AB5 ( R u } -RP ( 1 ) ) 

IFIRE. GT.ER.l6Xk)  EKmaXH*Kl 

SSQR«5SwR+R£»*2 

AE«AbS ( AZ ( 1 )-AP ( I ) ) 

I F ( AE.G) .ERMAXA  ) LRMAXA*AE 
U SSQA*SSWA*AE**2 

XMSWR*SSDR/FLOAT (N-J5J 
XMSWA*5b«A/FL0A I (N-lb  ) 

RMSWR*Sm(R  T I XMSUk  ) 

RMSWA*SwRT ( XMSUA  ) 

PRINT  b,  1N01SE, XMSQR,RMSUt\,ERMAXR, XMSUA, RmSUA  (ERMAXA 
10  CONTINUE 
STOP 
END 
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MAIN  PROGRAM  TSTINI 


DIMENSION  R!20U»  ,AZ!200)  ,xMR(200) ,XMA«200} 

DIMENSION  RPUOU)  ,RV(200)  ,RA(2QU>  , RES (230  I 
DIMENSION  AP(200).AV<2UQ),AA(200>*AESt2G0> 

REAL  K 
FORMAT  J ) 

format ( 5x * • g h k • .3F12.5) 

FORM AT ( 5X  » * R5TD  AS  T D ♦,  2F12.5) 

FORMAT ( 5X  » * INIT  *,15) 

FORMAT! /5X,  ' RESULTS  FOR  INOISE  » *,!$// 

* 20X.*  RANGE  MEAN  SQUARED  ERROR  * 'ifli.S/ 

* 2UX*»  RANGE  RMS  ERROR  « »,FI2.5/ 

* 20X,»  RANGE  MAXIMUM  ERROR  * *,FU«S/ 

* 20X » * ANGLE  MEAN  SQUARED  ERROR  » *,F12.S/ 

* xUXi*  ANGLE  RMS  ERROR  « «,F12*5/ 

* 20X,*  ANGLE  MAXIMUM  ERROR  ■ *,F12.5> 

FORMAT ( 5X  i ’ DEL  T »,F12#5) 

READ  t i DEL T 
PRINT  7 i DELT 
READ  1 i G • H • K 
PRINT  2 « G • H • K 
Rt.  AU  l,RSTU,A:>TO 
pH  i u T 3)RSTD»  A S T U 
READ  l.INIT 
PR  I NT  M , 1 N I T 

call  traj<r,az,n,deet> 

Oo  10  I NO  i St* i »3 
SSQR-U, 

SSUA*0» 

EnmAXr»0 • 

ERMAX A«0, 

CALL  NOISE(K,aZ*N»XNR,XMA*RSTDiASTD) 

CALL  GHKFIL! XMH  ,n  ,KP  ,KV  ,RA  ,DELT  .G  ,H  ,K  »KtS. I N I T » 

CALL  GHKFIL(XMA,N,AP,AV,AA|OELT,G»H,K,A£S»INIT) 

00  11  I ■ H . 1 5 
RE3 A05 ( R ( 1 )-RP«  I ) ) 

IF  ( RE.GT.ERMAXR  ) t.RnAXH«RE 

SSQR«SSQR+RE**2 

At3 A8S { A Z ( I ) -AP ( 1) ) 

IM  AE.GT.ERrtAXA)  ERtiAXA»AE 
SSQA«SSQA+At**2 
XMSQR«SSQR/FLOAT ( N- is  ) 

X W S Q A * S 5 Q A / F L 0 A T (N-iS! 

RMSyR3SQHT(XMSQK) 

RiiSuA3SuRT ( X M S 0 A ) 

PRINT  5»lN0lSt.XMSwR,RMSyt<,ERnAxR|XilSQA|RMS0AfENMAXA 
CONTINUE 
S I OP 
E NO 
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APPENDIX  B 

PROGRAMS  FOR  EVALUATION  OF  ADAPTIVE  FILTER 
This  appendix  contains  the  FORTRAN  listings  of  the  programs  which 
are  used  for  the  evaluation  of  the  adaptive  g-h-k  filter  presented  in 
Chapter  5.  These  programs  are  subroutine  GHKADA  and  main  program  TSTADA. 
Other  subroutines  required  for  this  study  are  TRAJ  and  NOISE  which  are 
listed  in  Appendix  B. 
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SUBROUTINE  GHKADA( XM,XP,XV.XA»THR) 

COMMON/FL T P AR/N, DEL T ,GN,HN,KN, WEIGHT , l STORE, GH,HH,XM, I TRANS 

01  HENS  I ON  XH(  l ) ,XP(  l ) ,XV (1)  ,XA( 1 ) 

DIMENSION  X5T(3,1Q) ,XS(3) 

DATA  XST /30*0 « / 

REAL  K * KN  * KM 

2 FORMAT ( 5X  * * MANEUVER  DETECTED  AT  IS  • **I5/) 

C INITIALIZATION 

XS(l)«(5**XM  ( 3 ! ♦ 2 « * X M (2)-XM  (1)1/6* 

XS ( 2 ) * ( XM  ( 3 ) “XM  ( 1) )/(2.*DELT) 

xsm«o* 

I S» 3 

XDET  *0* 

IrtAN*0 
G*GN 
H*HN 
K *KN 

DU  10  1*1 *3 
10  XbT(I,l)-XSm 
C PREDICTION  ALGORITHM 
IS  1S*ISM 

IFUS.GT.N)  RETURN 

XP( IS)*XS< 1 ) +D£LT*XSt  2 ) *XS( 3 ) • ( DEL Jt* 2/2*  ) 

X V ( IS)*XS(2)*DELT«XS(3) 

X A { l S ) *XS ( 3 ) 

RES  *XM  ( IS)-XP( IS) 

1F( IMAN.EQ*1.0R.IS*LT. IS)  GO  TO  12 
C MANEUVER  DETECTION 

XDET-WEIGHT*XDET*DELT*RES 
IF(ABS(XDET).GT*THR)  60  TO  11 

12  XSU  )*XP<  IS)«-G»Rf.S 

XS ( 2 ) *X V ( IS)*H*RES  /DELT 

XS ( 3 ) ■ X A ( 1S)+K*RES  ♦ 2 • / ( DELT • *2  ) 

C STORE  ISTORE  SMOOTHED  VALUES 
DO  13  1*1*3 

1 F ( 1 STORE « EG • l ) GO  TO  13 
DO  1R  J*I STORE  * 2 *•  l 
1 H XST  ( I » J ) *XST  ( 1 »J**1  ) 

13  XST ( 1 . 1 )*XS(  I ) 

IF( IMAN.EU.O)  GO  TO  15 
C CHECK  IF  TRANSIENT  PERlOO  IS  OVER 
IF( IS-IOET.LE. ITRANS ) GO  TO  15 
G • G N 
H*HN 
K*KN 
1MAN«0 
XUL  T *0  * 

GO  TO  15 

C REINITIALIZATION 
! I 10ET-1S 

PRINT  2,  IS 
I MAN* l 
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G»«M 

H*HM 

K*KM 

DU  16  1*1,3 

16  XSm*XST(  I , ISTORf  1 
00  17  1*1ST0R£.  1 }=•! 

PPX*XS<  1 l>0£4.T*XS(2)+XSC3i*CDELT**2/2*  1 
PtfXaXS(2!+0E|.T*XS(3) 

PAX*XS ( 3 1 

RESX  X *XM  UDET-I  M >-PPX 
XS< 1 )*PPX* G*RESXX 
XS(2)*PVX*H*RESXX/DELT 

1 7 XS13I»PAX*K*RESXX*2#/(0ELT**2) 

GO  TO  15 

END 


1 

2 

3 

H 


9 

21 

22 


MAIN  PROGRAM  TSTADA 

COMHON/FLTPaR/N»OE1.T»GN,HN,KN,WE1SHT,ISTORE»GM,HM,KMi1TR^NS 
Oi MENS  I ON  RP ( 2Q0 } ,RV ( 200 > , R* ( 2Q0 ) 

Dimension  ap(2Q0> ,av 12001 ,aai 200 > 

DIMENSION  XMR(200) .xMA(200) 

DIMENSION  R(20Q| ,aZS200) 

REAL  KM.KN 
FORMAT ( 1 
FORMAT  1 20X , * 

FORMAT  1 20X , * 

FORMAT  1 20X , * 

F0RMAT1SX,’ 


TRACKING  INTERVAL  • *,F12.S) 

COEFFICIENTS  1 NONMANEUVERING  • , 3 ( F 1 2 . S , 3X M 
COEFFICIENTS  1 MANEUVERING  • , 3 1 F 1 2 . 5 , 3X ) ) 
MAXIMUM  AND  RMS  ERRORS  •*/ 


2 0 X i 


• 2QX , 

• 2nx . * 


RANGE 

RANGE 

RANGE 

ANGLE 

ANGLE 

angle 


2QX, 

2QX,» 

• 20X,* 

• 20X,» 

FORMAT ( 2QX  , 

• 20X  i 

• 20X , 
FORMAT ( 20X , 

• 20  X , 

FORMAT  1 20X  , 

• 20X  , 

FORMAT { /5X  , 
FORMAT ( /5X  , 

F0HMATU5X, 

READ  1,0ELT 
PRINT  2 , DEL  T 
READ  l,GN,HN,KN 
PRINT  3 , GN  , HN  , K N 


■ »tFI2,5/ 

* *»F12.5/ 

* *,F12«5/ 

* * i F 1 2 • 5/ 

* »»FI2#5/ 

* * i F 12 • 5 ) 
WEIGHT  * 


MEAN  SQUARED  ERROR 
RMS  ERROR 
MAXIMUM  ERROR 
MEAN  SQUARED  ERROR 
RMS  ERROR 
MAXIMUM  ERROR 
EXPONENTIAL  f ILTER 
RANGE  THRESHOLD  » 

ANGLE  THRESHOLD  * 

SAMPLES  FOR  REINITIALIZATION 
SAMPLES  IN  TRANSIENT 
RANGE  STO,  DEVIATION  ■ 

ANGLE  STD,  DEVIATION  * 
RESULTS  FOR  INOISE  * * , 15/ ) 
RANGE  TRACKING  COMPLETE’/) 
ANGLE  TRACKING  COMPLETE*/) 


, F 1 2 , 5/ 
,F12.5/ 
t F 1 2 • 5 ) 
,15/ 
,IS> 

, F 1 2 » 5/ 

, F 1 2 « 5 ) 
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REAO  l,6M,HM,KH 
PRINT  9 , GM  * HM , KM 
REAO  I .WEIGHT  ,RTHR ,AThR 

PRINT  6, WEIGHT, RTHRiATHR 
REAO  I , I STORE  , I TRANS 
PRINT  7.IST0RE.ITRANS 
READ  1 »RSTD  » ASTD 
PRINT  8,RST0,ASTD 
CALt  TRAJ(R,AZ,N,DELT> 

DO  10  lNOISE-1,3 
PRINT  9 , 1 NO  1 SE 
SSqRbQ • 

5SQA-0, 

ERMAXR-O. 

ERMAX A»U . 

CALL  N01SUR,AZ,N,XMR»XMA,RSTD,ASTD) 

CALL  GHKA02(XMR,RP,RV,RA,RTHR) 

PRINT  21 

CALL  GHKA02(XHAiAP,AV,AA,ATHR) 

PRINT  22 

DO  11  I ■ 1 6 , N 

RE» ABS ( R ( I l-RPtll > 

IFIRE. GT.ERMAXR1  ERM AXR«RE 

SSQR“SSGR*RE**2 

AC- ABS ( AZ ( I )-AP( 1 ) 1 

I F { AE .GT .ErHAX A i ErMAXA-AE 

U SSQA»SSQA*AE**2 

XMSGR«SSGR/FL0AT(N-15> 

XMSQA-SSUA/FLOAT (N-15) 

RMSGR«S«RT ( XMSQR  > 

RMSGA-SWRT ( XMSGA ) 

PRINT  5,XmSuR,RMSQR,ERMAXR,XMSGA,RMSQa,ERMAXA 

10  CONTINUE 
STOP 
END 
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